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QUANTUM MONTE CARLO SIMULATION^ 

By Yazhen Wang 

University of Wisconsin-Madison 

Contemporary scientific studies often rely on tlie understanding 
of complex quantum systems via computer simulation. This paper 
initiates the statistical study of quantum simulation and proposes 
a Monte Carlo method for estimating analytically intractable quan- 
tities. We derive the bias and variance for the proposed Monte Carlo 
quantum simulation estimator and establish the asymptotic theory 
for the estimator. The theory is used to design a computational 
scheme for minimizing the mean square error of the estimator. 

1. Introduction. Computer-aided simulations of physical systems are wi- 
dely used in scientific and engineering studies such as aircraft and car design 
and nuclear explosion modeling. While the traditional simulation methods 
with the aid of classical computers based on transistors are to understand 
basic properties of materials, many contemporary simulations rely on un- 
derstanding quantum systems, such as those in bio-chemistry and nano- 
technology for the design of nano-materials and novel molecules. See Aspuru- 
Guzik et al. (2005), Kou (2009) and Waldner (2007). 

A quantum system is described by its state, which is often characterized 
by a vector in some complex Hilbert space. The number of complex num- 
bers required to characterize the quantum state normally grows exponen- 
tially with the size of the system, rather than linearly, as occurs in classical 
physical systems. Consequently, for a quantum system it takes an exponen- 
tial number of bits of memory on a classical computer to store its quantum 
state, and simulations of quantum systems via classic computers face great 
computational challenge. As quantum systems are able to store and keep 
track an exponential number of complex numbers and perform data ma- 
nipulations and calculations as the systems evolve, quantum computation 
and quantum information are to grapple with understanding how to take 
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advantage of the enormous information hidden in quantum systems and to 
harness the immense potential computational power of atoms and molecules 
for the purpose of information processing and computation. Quantum com- 
puters built upon quantum systems may excel in the simulation of naturally 
occurring quantum systems, where such quantum systems may be hard to 
simulate in an efficient manner by classical computers [Abrams and Lloyd 
(1997), Boghosian and Taylor (1998) and Zalka (1998)]. 

To the best of our knowledge, this paper is the first to introduce quantum 
computation and study quantum simulation in the statistical framework. 
Specifically, we will propose a Monte Carlo quantum simulation method 
for computing analytically intractable quantities and analyze approximation 
errors and random variations of the proposed Monte Carlo estimator. The 
theoretical analysis establishes a strategy to design an optimal scheme for 
utilizing computational resources in obtaining the Monte Carlo estimator. 

The rest of the paper proceeds as follows. Section 2 provides a brief review 
on quantum mechanics, quantum statistics and basic concepts of quantum 
computation. Section 3 proposes a Monte Carlo quantum simulation method 
and then presents the statistical analysis for the method. We derive the 
variance and bias for the proposed estimator and establish the strategy to 
allocate computational resources in the Monte Carlo quantum simulation for 
minimizing the mean square error of the estimator. A quantum simulation 
example is illustrated in Section 4. 

2. Brief background review. 

2.1. Quantum physics. Quantum mechanics describes phenomena at mi- 
croscopic level such as position and momentum of an individual particle like 
an atom or electron, spin of an electron, detection of light photons, and the 
emission and absorption of light by atoms. Unlike classical mechanics where 
measurements of quantities like position and momentum can be observed 
accurately, the quantum theory can only make statistical prediction about 
the results of the measurements performed. 

Mathematically quantum mechanics is usually described by a Hilbert 
space Ti and Hermitian (or self-adjoint) operators on Ti. As in quantum 
mechanics, we adopt standard Dirac notation |-), which is called a ket, to 
indicate that the object is an element inT-L. A quantum system is completely 
described by its state and the time evolution of the state. A state is often 
classified as a pure state or an ensemble of pure states that are easy to de- 
scribe by density operators. A pure state is a unit vector {i/j) in "H, which 
corresponds to a density operator p = \ip){tp\, the projection operator on 
An ensemble of pure states corresponds to the case that the quantum system 
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is in one of states IV'fc), k = 1, . . . , K , with probability pk being in state {ipk), 
and the corresponding density operator is 



K 

(1) P = ^Pk\if^k){il^k\- 

k=l 

Let {ipit)) be the state of the quantum system at time t. The states \^l'{tl)) 
and |V'(*2)) at ti and t2 are connected through \tp{t2)) = U{ti,t2)\'4'{ti)), 
where U{ti,t2) is a unitary operator depending only on time ti and t2. In 
fact, there exists a Hermitian operator H, which is known as the Hamilto- 
nian of the quantum system, such that U{ti,t2) = ex.p[—iH{t2 — ti)]. With 
Hamiltonian H, we may depict the continuous time evolution of |V'(i)) by 
Schrodinger's equation 

See Holevo (1982) and Sakurai (1995). 



2.2. Quantum probability. Quantum mechanics can be tested by check- 
ing its predications with experiments of performing measurements on quan- 
tum systems. The common quantum measurements are on observables such 
as position, momentum, spin and so on, where an observable is defined as 
a Hermitian operator on Hilbert space T-L. Consider an observable X with 
a discrete spectrum so that it can be written in a diagonal form 

p 

(3) X = y^XaQa, 

a=l 

where S K are eigenvalues of X and Qq are the corresponding one- 
dimensional projections onto the eigenvectors of X. Possible measurement 
outcomes of the observable are described by measure space {i},J-). For 
a quantum system with state p, the result of the measurement is random 
with probability distribution Pp over ($7,7^). We denote by X the result 
of the measurement of observable X given by (3). The result X is a ran- 
dom variable and takes values in = {xi,X2, ■ ■ ■}■ With a quantum sys- 
tem prepared in the state p, the result X has a probability distribution 
Pp[X = Xa] = tr(/3Qa). With the probability distribution Pp, we can easily 
derive the expectation and variance of X in the state 

p 

Ep[X] = tr(pX) = Y,^aPp[X = Xa] = Ep^iX), 

a=l 

Varp[X]=tr[pX2]-[tr(pX)]2. 
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Measuring the outcomes of observable X will alter the state of the quan- 
tum system. If the state of the quantum system is p immediately before the 
measurement, then the probability that the result Xa occurs is Pp[X = Xa\ = 
tr(/>Qa) and the state of the system after the measurement result Xa is equal 
to Qa/oQa/ tr(Qa/3Qa). Similarly, using the spectral theory of self-adjoint 
operators, we may describe observables with continuous spectrum and con- 
tinuous measurement outcomes. See Barndorff-Nielsen, Gill and Jupp (2003) 
and Holevo (1982). 

2.3. Quantum computation. Quantum systems can be simulated via com- 
puters, but quantum simulation requires enormous computational resources. 
Classic computers may have great difficulty to efficiently simulate general 
quantum systems, while quantum computers built upon quantum systems 
are ideal for quantum simulation. 

Analog to the fundamental concept of the bit in classical computation 
and classical information, we have quantum bit in quantum computation 
and quantum information and call it qubit for short. Just like a classical 
bit with state either or 1, a qubit has states |0) and However, there 
is a real difference between a bit and a qubit. Besides states |0) and 
a qubit can also take states as their superpositions, which are the linear 
combinations of |0) and 

|^) = ao|0)+Qi|l), 

where complex numbers ao and ai are called amplitudes satisfying |aoP + 
|qiP = 1. In other words, the states of a qubit are unit vectors in a two- 
dimensional complex vector space, and states |0) and |1) consist of an or- 
thonormal basis for the space and are often referred to as computational 
basis states. The qubit is the simplest quantum system. Unlike a classical 
bit which can be examined to determine whether it is in the state or 1, 
for a qubit we can not determine its state and find the values of oq ^-nd ai 
by examining it. Quantum mechanics shows that we can measure a qubit 
and obtain either the result 0, with probability |aoP, or the result 1, with 
probability |aip. A qubit can be actually realized as physical objects in 
many different physical systems, such as the two different polarizations of 
a photon, the alignment of a nuclear spin in a uniform magnetic field or 
two states of an electron orbiting a single atom. In the atom model case, we 
may correspond |0) and |1) with the so-called "ground" or "excited" states 
of the electron, respectively. As the atom is shined by light with suitable 
energy and for a proper amount of time, we can move the electron from 
the |0) state to the |1) state and vice versa. Moreover, by shortening the 
length of time shining the light on the atom, we can move an electron ini- 
tially in the state |0) to "halfway" between |0) and say, into a state 
|+) = (|0) + |l))/^/2. 
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Like classical bits, we may consider multiple qubits. The states of two 
qubits are unit vectors in a four-dimensional complex vector space, with 
four computational basis states labeled by |00), |01), |10) and |11). In gen- 
eral, a system of b qubits has 2'' computational basis states of the form 
\x1X2 • • • Xb), Xj = or 1, j = 1, . . . , b, that generate a 2^-dimensional com- 
plex vector space, and a superposition state in the system is specified by 2* 
amplitudes. As 2'^ increases exponentially in b, it is easy for such a system 
to have an enormously big vector space. A quantum system consisting of 
even a few dozens of "qubits" will strain the resources of even the largest 
supercomputers. Consider a system with 50 qubits. 2^^^ ~ 10^^ complex am- 
plitudes are needed to depict its quantum state. With 128 bits of precision, 
it requires approximately 32 thousand terabytes of information to store all 
10^^ complex amplitudes. Had Moore's law continued on schedule, such stor- 
age capacity would be available in supercomputers during the second decade 
of the twenty-first century. A system with b = 500 qubits has the number of 
amplitudes larger than the estimated number of atoms in the universal. It is 
unimaginable to store all 2^'''^ complex numbers in any classical computers. 
In principle, a quantum system with only a few hundred atoms can manage 
such an enormous amount of data and execute calculations as the system 
evolves. Quantum computation and quantum information are to find ways 
to utilize the immense potential computational power in quantum systems. 
See Clarke and Wilhelm (2008), Deutsch (1985), DiCarlo et al. (2009), Feyn- 
man (1982), Lloyd (1996), DiVincenzo (1995), Nielsen and Chuang (2000) 
and Shor (1994). 

3. Statistical analysis of quantum simulation. The key for the simula- 
tion of a quantum system lies in the solution of Schrodinger's equation (2) 
which governs the dynamic evolution of the system, and the quantum sim- 
ulation can be done via either classic computing or quantum computing. 
Schrodinger's equation for a typical Hamiltonian with real particles usually 
consists of elliptical differential equations, where each differential equation 
can be easily simulated by a classical computer. The real challenge in stim- 
ulating a quantum system is to solve the exponential number of such differ- 
ential equations. Consider a quantum system that is described by b qubits. 
As b qubits have 2'' amplitudes, for stimulating the dynamic behavior of 5 
qubits evolving according to Schrodinger's equation, a system of 2^ differ- 
ential equations must be solved. Because of the exponential growth in the 
number of differential equations, simulating quantum systems by classical 
computers is feasible only for special cases where insightful approximations 
are available to dramatically reduce the effective number of differential equa- 
tions involved. Quantum computers may be ideal for the simulation of natu- 
rally occurring quantum systems. See Abrams and Lloyd (1997), Boghosian 
and Taylor (1998), Feynman (1982), Lloyd (1996) and Zalka (1998). Whether 
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quantum simulation is via classic computing or quantum computing, its sta- 
tistical aspect essentially remains the same. 

3.1. Simulate a quantum system. The heart of quantum simulation is to 
solve Schrodinger equation (2) which has the solution 

(4) |^(t)) = e-^*|^(to)>. 

Numerical evaluation of e~*^* is needed. The Hamiltonian H is usually ex- 
ponentially large and extremely difficult to exponentiate. The common ap- 
proach in numerical analysis that uses the first-order linear approximation, 
1 — iH6, of e~*^(*+'') — e~*^* often yields unsatisfactory numerical solutions. 

Efficient evaluation of the solutions (4) with high order approximation 
exists for many classes of Hamiltonians. For most physical systems the 
Hamiltonians involve only location interactions, which originate from the 
fact that most interactions fall off with increasing distance or difference in 
energy. Specifically, a system of a particles in a d-dimensional space often 
has a Hamiltonian of the form 

L 

(5) H = Y,Hi, 

e=i 

where L is a polynomial in a -|- d, and each Hi acts on a small subsystem of 
finite size free from a and d. Typical examples of the terms Hi are one-body 
Hamiltonians and two-body interactions such as the ones in the Hubbard and 
Ising models [Altland and Simons (2006) and Dziarmaga (2005)]. As a re- 
sult, e~^^'^ is easy to compute numerically, although e"*^*^ is very hard to 
evaluate. Because Hi and H^ are noncommutable, e"*^"^ / g-^H^S . . . ^-iH^S ^ 
Using the Trotter formula [Kato (1978), Sornborger and Stewart (1999) and 
Trotter (1959)], we approximate e~^^^ by Us which requires only the eval- 
uation of each e~*^^^, where 

Suppose that the quantum system starts at to with initial state \^jJ{to)) 
and ends at final time T. For an integer m, let 5 = T/m and tj = j6, j = 
0,1, ... ,m. The quantum simulation is to apply approximation Us of e~*^^ to 
solutions (4) at tj iteratively and generate approximate solutions for |^(tj)). 
Denote by \^Jj{tj)) the state at tj obtained from the quantum simulation as 
an approximation of the true state \ip{tj)) at tj. Then for j = 1, . . . ,m, 

m,)) =e-'^'mj.,)) =e-'^^'mo)), 

(7) , , . 

m,)) = usmt,^i)) = ujmto))- 
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If the initial state of the quantum simulation is a pure state \'tlj{tQ)), then 
the true final state and the simulated final state are also pure states \tp{tm)) 
and \ip{tm)), respectively, with corresponding density operators 

p{to) = mtm)){^{trn)\, /5 = | . 

When the initial state is an ensemble of pure states, with probability pk 
being pure state {ipkito)), k = 1, . . . , K, and corresponding density operator 

K 

p{to) = ^Pk\'4'k{to)){ipkito)\, 
k=l 

then at time tj the true state and the simulated state are also ensembles of 
pure states with respective density operators 

K 

P(h) = ^Pk\'4'k{tj)){i^kitj)\, 
k=l 
K 

P(h) = ^Pk\i'k{tj)){i>kitj)\, 
k=l 

where for fc = 1, . . . , K, 

\Mtj)) = U5\Mtj-i)) = ui\Mto))- 

3.2. Monte Carlo quantum simulation. Quantum simulation provides an 
excellent way for the study of complex phenomena in physical and biology 
systems and evaluating hard-to-obtain quantities in the system. Examples 
include the dielectric constant, the mass of the proton, conductivity, mag- 
netic susceptibility of materials and molecules in biological systems. As the 
results of quantum measurement outcome are random, repeated measure- 
ments need to be performed in order to obtain reliable estimators of the 
quantities. In the quantum setup, a quantity of interest is of the form 

(8) e = E,iX)=Tr{Xp) = Ep^iX), 

where X is an observable, X is its measurement result, and p is the state of 
the quantum system under which we perform the measurements and evaluate 
the quantity 9. 

A Monte Carlo quantum simulation method is designed to estimate 9 as 
follows. We prepare the quantum system at initial state p{to) and make it 
evolve to final state p{tm) = P- The quantum simulation procedure described 
in Section 3.1 is used to simulate the evolutions of the quantum system from 
initial state p{to) to final state p{tm) according to Schrodinger's equation (4) 
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with some Hamiltonian H of the form given by (5). We repeatedly perform 
the measurements of observable X in such n identically simulated quan- 
tum systems at the simulated final state and obtain measurement results 
Xi,. . . , Xn- We estimate 9 defined in (8) by 

(9) ^=lt^''^- 

Since measurements and state approximations in quantum simulation are 
involved with random fluctuations and systematic errors, as a Monte Carlo 
estimator of 9 has variance and bias. We use mean square error (MSE) crite- 
rion to gauge its performance. The Monte Carlo quantum simulation method 
for obtaining estimator 9 requires to repeat the simulation of the quantum 
system n times, and each simulation needs to calculate m approximations of 
states at tj, j = 1, . . . ,m. The whole Monte Carlo quantum simulation pro- 
cedure needs to make total N = mn state approximations for the quantum 
system. One important problem is to determine the strategy to allocate m 
and n with given N = mn so that the MSE of 9 is minimized. We derive the 
MSE of 9 and study the problem in the following theorem. 



Theorem 3.1. For a quantum system evolving in time interval [0, T], 
assume that its Hamiltonian H and observable X satisfy (5) in Section 3.1 
and (11) in Section 3.3, respectively. Then 

m-9f]<^+^ = ^^ + ''''" 



n m 



4 r4 ' 



where Ci and C2 are generic constants free from m and n. Thus, when 
n = C\vr& IC2 and m = {C2/CiY^^N^/^ , the MSE bound is asymptotically 
minimized and 

E[{9-9f]<Ct"'Cl'''N-^l\ 



The theorem indicates that the MSE of 9 is of order C\n~^ + C2m~'^, 
where from the proof of the theorem in Section 3.3 below we see that Ci is 
the variance of X and C2 is the difference of the expectations of X under 
true state p{tm) and the simulated state p{tm)- As the Monte Carlo quan- 
tum simulation procedure performs n repeated simulations of the quantum 
system with m state approximations for each simulation, if we have the 
computational capacity of carrying out a total of N = mn state approxima- 
tions in the Monte Carlo quantum simulation, the theorem provides optimal 
strategy for the allocation of m and n that minimizes the MSE of 9. 
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3.3. Proof of Theorem 3.1. As usual, the MSE analysis involves deriving 
its variance and bias, 

(10) E[{B-ef]=\ixv{e) + {Ee-ef. 

We need to fix some notation to facilitate further analysis. The target 9 
is defined under the true final state p{tm), while the quantum simulation 
is under approximate final state p{tm) which is close to p{tm)- To make 
the problem realistic, we impose the following assumption to ensure that 
observable X behaves well in states close to the true final state. With the 
true final state of the form 

K 

P{tm) = y^Pfc|'0fc(^m.))(V^fc(^m)|, 
k=l 

we assume that for some small ?y > 

(11) max sup{||X|0)||,|||0)-|V^fc(t^))|| <77}<cx). 

l<k<K 

The condition is to ensure that observable X has two finite moments under 
states in a small neighborhood of the true state p{tm) of the quantum system. 

Since a simple conditional argument will reduce the proof from the gen- 
eral ensemble state case to the pure state case, for simplicity, we consider 
the Monte Carlo study with pure states. The state used in (8) is the pure 
state \tp{tm)) or p= \^p{tm)){'>P{tm)\, under which observable X and its mea- 
surement result X are considered. The measurement results Xi, . . . ,Xn are 
obtained from the quantum simulation under the simulated state |V'(tm.)) 
OT p= \ip{tm)){'4'itm)\- Therefore, to analyze the bias and variance, we need 
to evaluate the expectations and variances of Xi under p but compute the 
corresponding quantities of X under p. 

The bias E'^ — = Tr(Xyo) — Tr(Xp) is due to the differences between \ijj{tj)) 
obtained in the quantum simulation and the true quantum states \'tp{tj)). 
We will prove in Proposition 3.1 below 

(12) \mu)-mtm))\\<c6\ 

Thus, we derive the bias 

\Ee -e\ = \ Tr(Xp) - Tr(Xp)| 

= \{^P{t„^)m{tra)) " (V'(U) |X| V'(trn)) | 

(13) <KV^(tm)-V'(tm)|X|^(i^))| + |(^(i^)|X|V;(tm)-V'(*m))| 

<iiiv;(t™)-V'(t„^))ii(iixiv;(tm)>ii + 11X1^^(^)11) 
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where the last inequahty is due to (12) and condition (11). The variance of 6 
is easy to obtain 

(14) Var(0) = - Var(Xi) = -{Tr(X'^p) - [Tr(Xp)]^}. 

n n 

As we have shown above, the Tr(Xp) approach to 6 = Tr(Xp) as m — ?■ cxd 
or, equivalently, 5— )-0. As for Tr(X^/)), 

I Tr{X.^p)\ = \{i,it^)\X.^mU)\ = 

whose finiteness is a consequence of (12) and condition (11). Collecting to- 
gether (10), (13) and (14), we conclude 

n n 

which is asymptotically minimized when n ~ ~ N^^^ . To complete the 
proof of Theorem 3.1, we show (12) in the rest of the section. 

The quantum simulation uses \tp(tj)) to approximate the true quantum 
states 11^(1 j)) and thus results in approximation errors. We define the follow- 
ing quantity to measure the approximation errors. Suppose Ui and U2 are 
two unitary transformations, the operator norm of the difference between Ui 
and U2, 

r([/i,C/2)= max||([/i-[/2)|</.)||, 

is used to measure the closeness of C/i and U2- 

In the quantum simulation scheme, we approximate e~^^^^ by Ug defined 
in (6). Naturally we use T{U^ ,e~^^^^) to gauge the approximation errors 
in the quantum simulation. Below we derive the order in terms of 5 for 
approximation errors — e~^^^^ in the quantum simulation. 

Proposition 3.1. The following inequality holds uniformly for j = 
l,...,m: 

T{uie-'^^^)<CL5^. 

Proof. First we prove the inequality for j = 1- For the case of L = 2, 
H = Hi + H2 ■ Expanding exponential functions of Hi and using simple al- 
gebraic manipulation, we have 

2 8 
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e"'^'''^''^^^ = 1- i{Hi + H2)6 - ^(^1 + H2fS^ + 0{6^) 

= 1- i{Hi + H2)5 - ^6\HiH2 + H2H1 + Hf + Hi) + ©(d^) 
= (l- \h,5 - \h15^^ [i - ^H26 - 
x[i-Ih,6-\h15'Yo{5') 

For general L, let H* = Yh=^ ^i- Then H = HI and H* = Hj + H*^^ for 
j = 1, . . . ,L — 1. We repeatedly apply the above result for the case of L = 2 
to the case of H* = Hj + H*^-^ and obtain 

^ ^-iHiS/2^-iH^5^-iHiS/2 _^ q^^3^ 

^ ^-iHiS/2^-i{H2+H^,)5^-iHi5/2 _^ q^j3^ 

^ ^-iHi5/2^-iH25/2^-iH;5^-iH25/2^-iHi5/2 _^ Q(2(5^) 

= ... = C/5 + 0(L53), 

which implies the inequality for j = 1. 

Second we show the inequality for j = 2. Let Vs = e~*^^. For any state |(/)), 

\\{u!-vi)m<\m!-usvs)m\ + \mvs-vi)m\ 

< \\Us{Us - Vs)\4>)\\ + \\{Us - Vs)Vs\4>)\\ 

<\\iUs-vs)m\ + ms-vsM')\\, 

where {(j)') = Vs\(j)). Since Vs is unitary, = |||0)||- On both sides of the 

above inequality we take the maximum over all (p with \\\(j))\\ = 1 and obtain 

r(C/|,F/)= max \\(ui-Vs'm\\ 

lll<p)ll=i 

<2ma^\\{Us-Vs)m\ 

= 2T{Us,Vs)<CL5^ 

where the last inequality is from the proved case of j = 1. 

The result for general j follows by induction. Since Us and Vs are unitary, 
we repeatedly apply the same technique for proving the case of j = 2 as 
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follows. For |(/>) with |||0)|| = 1, 



wi - vimw < \\{ui - uiWs)m + m^-Ws - vi)m\ 
< \\ut\us - vsMw + \\{ui~' - vr'wmi 
<\\{us-vsm\ + \\iur'-vr'M')i 



where {(j)') =Vs\(j)). Taking the maximum over all \(j)) with \\\(j))\\ = 1 on both 
sides of the above inequality, we get 



where the last inequality is from m5 = T and the proved case of j = 2. □ 

4. An example. There are a few interesting and realistic quantum sys- 
tems such as the quantum Ising model and simple harmonic oscillator for 
which some analytic solutions are available [Dziarmaga (2005) and Sakurai 
(1995)]. In this section we illustrate the Monte Carlo quantum simulation 
method with simple harmonic oscillator. 

4.1. Three-dimensional isotropic harmonic oscillator. We consider a quan- 
tum system oi d/3 particles with three-dimensional isotropic harmonic oscil- 
lator. This is a d-dimensional quantum system. By using the natural scales of 
length and energy in terms of particle mass, angular frequency and Planck's 
constant, we have the following simple expression for the Hamiltonian of the 
system: 



and for £ = 1, . . . ,d/3, {x3e-2, x^£_i, x^i) specify the position coordinates of 
the ith. particle in M^. Hamiltonian H can be written as a sum of d local 



and have expression Hj = A +1/2, where and A are creation and 
annihilation operators given below. 



mivi) < r{Us,Vs) + r{Wf\ y/-^) < • • • < jr{Us, Vs) 

<mr{Us,Vs)<CL6'^, 




d d 





A^ = T Vj)/V2, [A-,A+] = A-A+- A+A- = I. 
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As shown in Dziarmaga (2005), the Hamiltonian of the quantum Ising model 
can also be expressed by similar product of creation and annihilation oper- 
ators. 

Operator A~^A~ has eigenvalues k for A: = 0,1,..., with eigenfunctions 
defined by normalized Hermite polynomials, 

In fact, it can be directly verified that [x'^hk{x) — h'^{x)]/2 = {k + l/2)hk{x) 
and, thus, hk{x) are eigenfunctions of Hj corresponding to eigenvalues k-\- 
1/2. As Hamiltonian H \s a. sum of d one-dimensional harmonic oscillators, 
eigenfunctions of Hamiltonian H are given by 

d 

^fcW = n^'^j(^j)' k = {ki,...,kd)\kj =0,1,..., 
i=i 

with corresponding eigenvalues Yl'j=i % + ^Z^- 

4.2. Quantum simulation. To make the quantum simulation manageable 
computationally, we consider the simulation of the following six-dimensional 
quantum system described by 12 qubits. It requires a Hilbert space of di- 
mension 2^^ = 4096 to accommodate the quantum system. We use the first 
four eigenfunctions of a one-dimensional harmonic oscillator to form two 
qubit states in each dimension. With the six sets of the four eigenfunctions, 
we obtain 4096 eigenfunctions of product form and generate a Hilbert space 
of dimension 4096 to accommodate the 12 qubit quantum system. To define 
and code the 12 qubits through the eigenfunctions, let z = (zi, . . . , zi2)^ with 
Zj = or 1, and k = {ki, . . . , k^)^ = {zi,. . . , zq)^ + 2{z7, . . . , zu)^ where f de- 
notes the transpose of a vector. The coordinates of k take four integer values 
from to 3. We identify qubit state \zi - ■ ■ Z12) with eigenfunction /i^(x). The 
quantum system is governed by Hamiltonian H and evolves in time inter- 
val [0, 1]. Let Vs = e~^^^. We illustrate the quantum simulation by approxi- 
mating Vs with 



Assume that the quantum system has an initial state at to = 0: 

1 ^ 3 



1 ^ 1 ^ 
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and final true state \'^m) at = 1) where for j = 1, . . . ,m, 

= e-^*^' Iv^o) = V,^M = ^ E e-(^^+-+^«+2..+...+2...+3)*, . . . ^^^^ 

kj=0 

The approximation states in the quantum simulation are 

Consider a path-dependent observable 
1 1 

z,=0 

where Q is a projection operator, z = (zi, . . . , 2:12)^, z{b) = Yl]=i ■^j^"'"^ 
responds to the number with binary representation zi - ■ ■ Z12, and 

g-«-ffiz(b) ^2,) = e"*^^!"' |-26+227 + ---+22l2+3)t^(6) |^_^ . . . 

We compute tr(Xp) and tr(X^p) as follows: 

1 

2^12) 



1 1 

Xlv^m) = ^ E h + 2^:7 H 
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1 ^ 

= ^ E + • • • + ^6 + 3)e-*('=^+-+'=«+3)/ifc^(xi) • • • /ifc,(x6), 

kj=0 

1 ^ 

= tr(Xp) = (<^^|X|(^„) = ^^^^ E + • ■ • + ^6 + 3) = 0.6, 

1 ^ 

tr(XV) = ((/^mlX^lv^™) = E + • • • + ^6 + 3)2 = 0.37875, 

Var(^) = i[tr(XV)-e2] = ^^^. 
n n 

Hence, we obtain the following expression for the MSE of 9: 

01 87^ 

MSE = + ((^^|X|(^^^) - 0.6)2. 

n 
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We need to numerically compute {ipm\'^\'^m) for the MSE evaluation. As Hi 
and Hj are commutable, e~*^-'^ = e~^^^^^ ■ ■ ■ e~^^'^^'^, and 

6 

The numerical method in Zalka (1998) can be used to evaluate e~^^'^^hk^ (xi) 
by repeatedly applying 

£=1 

We approximate {{(pm\^\'f>m) — 0.6)^ for = 5000 and 5 ranging from to 
0.01 and then evaluate MSE. The resulting MSE as a function of 5 decreases 
for 5 from to 0.0035 and then starts to increase. Its unique minimum is 
achieved at 5 = 0.0035, which corresponds to m = 277 and n = 18. Thus, 
with total 5000 times of state approximations allowed in the Monte Carlo 
quantum simulation for estimating 0, the Monte Carlo strategy to minimize 
the MSE of 9 is to take 5 = 0.0035 in the quantum simulation scheme and 
repeatedly simulate the quantum system 18 times. 

Acknowledgments. The author thanks editor Samuel Kou, the Associate 
Editor and an anonymous referee for helpful comments and suggestions. 
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